Time-dependent variational principle for quantum lattices 
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We develop a new algorithm based on the time-dependent variational principle applied to matrix prod- 
uct states to efficiently simulate the real- and imaginary time dynamics for infinite one-dimensional quan- 
tum lattice systems. This procedure: (1) is argued to be optimal; (2) does not rely on the Trotter decom- 
position and thus has no Trotter error; (3) explicitly preserves all symmetries and conservation laws; and 
(4) has low computational complexity. The algorithm is illustrated using both an imaginary time and a 
real-time example. 



The density-matrix renormalization group (DMRG) is 
arguably the most powerful tool available for the study of 
one-dimensional strongly interacting quantum lattice sys- 
tems [1]. The DMRG — now understood as an application 
of the variational principle to matrix product states (MPS) 
[2] — was originally conceived as a method to calculate 
ground-state properties. However, there has been a recent 
explosion of activity, spurred by insights from quantum in- 
formation theory, in developing powerful extensions allow- 
ing the study of, e.g., finite-temperature properties, higher- 
dimensional systems, and nonequilibrium physics via real- 
time evolution [3]. The simulation of nonequilibrium prop- 
erties with the DMRG was first attempted in [4], but mod- 
ern implementations are based on the time-evolving block 
decimation algorithm (TEBD) and relatives [5]. 

At the core of a TEBD algorithm lies the Lie-Trotter de- 
composition for the propagator exp(idtH), which splits it 
into a product of local unitaries. This product can then be 
dealt with in a parallelised and efficient way: when applied 
to an MPS one obtains another MPS with larger bond di- 
mension. To proceed one then truncates the MPS descrip- 
tion by discarding irrelevant variational parameters. This 
is such a flexible idea that it has allowed even the study 
of the dynamics of infinite translation-invariant lattice sys- 
tems via the iTEBD [6]. Despite its success the TEBD 
has some drawbacks: (1) the truncation step may not be 
optimal; (2) conservation laws, e.g. energy conservation, 
may be broken; and (3) symmetries, e.g., translation invari- 
ance, are broken (although translation invariance by two- 
site shifts is retained for nearest neighbor Hamiltonians). 
The problem is that when the Lie-Trotter step is applied to 
the state — stored as an MPS — it leaves the variational 
manifold and a representative from the manifold must be 
found that best approximates the new time-evolved state. 
There are a variety of ways to do this based on diverse 
distance measures for quantum states but implementations 
become awkward when symmetries and conservation laws 
are brought into account. 

In this Letter we introduce a new algorithm to solve the 
aforementioned problems — intrinsic to the the TEBD — 



without an appreciable increase in computational cost. The 
resulting imaginary time algorithm quickly converges to- 
wards the globally best uniform MPS (uMPS) approxima- 
tion for translational invariant ground states of strongly 
correlated lattice Hamiltonians, and the corresponding 
real-time evolution evolves an initial state without violat- 
ing energy conservation for constant Hamiltonians, or the 
conservation of any other quantities dictated by symmetry. 
The complexity of our approach can be made to scale as 
D 3 , comparable with current implementations, where D is 
the bond dimension of the uMPS. 

We now introduce the variational manifold .M u mps of 
uniform MPS for an infinite lattice of spin-d/2 degrees of 
freedom, parameterized via 

m^)) =T l d {sk} =A(UnezA s -)v R \s) , (1) 

where |s) = |. . . SiS 2 ■ ■ • ) and v L and v R are two D- 
dimensional vectors, which are presently argued to be ir- 
relevant. The variational parameters A comprise the set of 
D x D matrices A s (s = 1, 2, . . . , d) and are denoted via 
a dD 2 vector with entries A 1 = A s a ~, with i = (a, s, f3) 
a collective index. The uMPS variational manifold has a 
gauge invariance: replacing A s i— >■ GA S G~ 1 for invert - 
ible G results in an identical state. We do not fix the gauge 
and simply assume that A s are completely general com- 
plex matrices. We do, however, assume that the transfer 
matrix E = 2~2t=i A 8 ® A s nas precisely one eigenvalue 1 
with corresponding left and right eigenvectors (l\ and |r) of 
length D 2 , to which we can associate DxD matrices I and 
r, respectively, by simply reshaping them. These matrices 
are Hermitian and positive and assumed to have full rank. 
We choose the normalization so that (l\r) = Tr(Zr) = 1. 
In addition, we assume that all other eigenvalues of E lie 
strictly within the unit circle, i.e. the spectral radius of 
E — \r)(l\ is smaller than 1. These conditions allow one to 
write for any local operator O acting on n contiguous sites: 

0(A,A) = W(A)\d\iP(A)) / (iP(A)\iP(A)) = 
(l\ £ 0*,., n , s „ ,M S1 ■ ■ ■ A s ") ® (A* 1 • • • A*") \r). 

s,t=l 
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The boundary vectors v L and v r do not feature in normal- 
ized expectation values and thus do not contain any varia- 
tional degrees of freedom. 

Denote a translation invariant nearest-neighbour Hamil- 
tonian as H = En£Z T n h n T~ n , where T is the shift op- 
erator and h acts non-trivially only on sites zero and one. 
We now try to approximate the time evolution generated by 
H of a uMPS \ip(A)) without ever leaving the variational 
manifold of uMPS with fixed bond dimension D, by in- 
troducing a time-dependent parameterisation A{t). Inser- 
tion into the time-dependent Schrodinger equation results 
in A 1 \d t ^{A(t))} = -iH \ip(A(t))), where we denote d, 
for d/dA\ Whereas the left hand side (LHS) is a linear 
combination of the tangent vectors \diifj(A(t))) that span 
the tangent plane T A Ai uM ps, the right hand side (RHS) is 
a general vector in Hilbert space and this equation does not 
have an exact solution for A 1 . The best approximation is 
obtained by minimizing 

||i i |^W))) + i^lV'W)))ll- 

The minimum is found by orthogonally projecting the evo- 
lution vector H \ip(A(t))) onto the tangent plane, as illus- 
trated in Fig. 1 . The resulting solution is determined by 



(c^/W> A 1 = -i{d?l>\HW), 



where the argument A(t) in every vector has been omit- 
ted for the sake of brevity. The LHS of Eq. (2) 
contains the dD 2 x dD 2 Gram matrix of the tan- 



gent vectors G^j (A, A) 



(dMA)\d 3 ij(A)). Ex- 



pressions for this Gram matrix and the vector in 
the RHS of Eq. (2) are best derived using the ex- 
plicit form for the tangent vector B % \ diip{A)) = 

E„ e z T n E? Sfc }=i A>->B>°A» ■ ■ ■ ) v R |s), and 

are given by 



B< G- lJ B j 



|Z| 



(l\E*\r) 



+ - E)-'E B A \r) + {l\E B A {\ - E)^E A \ 

+ (|Z - l)\{l\Ei\r){l\E B A \r) 

I?(m&W) = \n[{l\Hi A \r) + {l\H£i\r) 
+ {l\Hft{l - E)^E A \r) + (l\E A (l - E)-'HH\r) 



where E A 



+ - 2){l\E A \r){l\H AA \r) 



E«=i A s (g>B (note the identity E = E A ) 



and H AB = Ei*, u ,„=i (s,t\h\u,v) (A^B") ® (C D ). 
In these expressions, (1 — E)^ 1 should be interpreted as 
the pseudo-inverse of (1 — E), i.e. it produces zero when 
acting on the left or right eigenvector of E with eigenvalue 
1: (Z|(l - E)- 1 = = (1 - E)- X \r). The overall fac- 
tors |Z| are a consequence of the infinite volume of our 
system and cancel, as they appear both in the LHS and 




d 2 i>(A)) 



FIG. 1. An illustration of our construction: the wireframe sur- 
face represents the variational manifold M = A1 u mps embedded 
in state space, with the black dot a point representing a uMPS 
\ip{A)). The rotated gray square represents the tangent plane 
T A M to M in \ip(A)), with two generally non-orthogonal co- 
ordinate axes \d\%p{A)) and |c>2(^4)} displayed as dotted lines. 
The arrow with solid head is the direction iH \ip(A)} of time evo- 
lution, and the arrow with open head represents the vector that 
best approximates iH \ip{A)) within the tangent plane. The gray 
curve is the optimal path \ip(A(t))) which follows the vector field 
generated by these vectors with open head throughout M. 



RHS of Eq. (2). The additional divergent terms on the last 
line of the brackets would disappear if we restricted our- 
selves to tangent vectors that are orthogonal to the uMPS 

= 0. 



(2 ) itself, such that (ip(A)\diip(A)} B l = |Z|(i|£7j 



Indeed, the tangent plane contains the state itself, since 
A i \di^{A)) = |Z| [ip{A)). However, a change in that di- 
rection would change the norm or phase of the state, which 
is not a desired effect. 

This construction can also be derived from an action 
principle and is known as the time-dependent variational 
principle (TDVP) [7, 8]. The resulting TDVP equations 
[Eq. (2)] can be shown to be sympletic [9]. Hence they 
respect energy conservation as well as conservation of all 
constants of motion, such as the expectation value of gen- 
erators of symmetries. Since only expectation values occur 
in the equations of motion, one can use techniques familiar 
from DMRG, including the decomposition of the matrices 
A 1 into irreducible representations of the relevant symme- 
try group. Further, this approach is manifestly translation 
invariant. For time -reversal invariant operators the TDVP 
equations are also invariant under time reversal (see [10] 
for a Trotter-based approach that recovers time reversal in- 
variance). This approach does not require any truncation 
and is thus globally optimal within the manifold A^ u mps- 

Constructing the relevant quantities and solving Eq. (2) 
for A 1 involve operations with a computational complexity 
of 0(D e ). Using an iterative method to implement (1 — 
E)^ 1 and then solving for A 1 can reduce this to 0(D 3 ). 
However, the matrix G^j is not invertible: because of the 
gauge invariance in the (u)MPS parameterisation, not all 
dD 2 tangent vectors are linearly independent. Defining 
the action of a 1 -parameter group of gauge transformations 
G(e) = exp(eX) as A s (e) = G(e)A s G(e)-\ we ob- 
tain that dA s jde = XA S — A S X. Because of gauge in- 
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variance, there is no corresponding change in \^(A(e))) 
and thus d\ip(A(e))} /de = (dA l / 'de) \d^) = 0. In- 
deed, any vector B l x defined by B s x = XA S — A S X pro- 
duces a zero norm state, evident when introducing it into 
the explicit form of B % \diip(A)). The vectors B x thus 
span the null space of Gjj. Any vector B in the tangent 
plane is gauge equivalent to B' = B + B x , VX G C DxD . 
There are D 2 — 1 linearly independent choices of B x , as 
we can easily prove by noting that B x = requires that 

ELi( AS V lB x = = ^ d a=1 (A')UXA'-lX. Since E 
has a single eigenvalue 1, and I has full rank, the only solu- 
tion to this equation is X = 1. In order to invert G^j, we 
fix the gauge which eliminates D 2 — 1 components of B. 
Norm preservation (i.e. (l\E^\r)) fixes one more compo- 
nent, resulting in a (d — 1)D 2 dimensional tangent plane. 

While there are a variety of ways to fix the gauge of vec- 
tors in the tangent plane, different choices result in different 
effective Gram matrices with different condition numbers. 
By using the gauge fixing condition (l\E% = — which 
also includes norm preservation and imposes the condition 
that the eigenvalue and left eigenvector of the transfer ma- 
trix do not change to first order — the effective Gram ma- 
trix reduces to B Jl G hj B j = \Z\(l\E§,\r) and all non- 
local contributions are thus effectively canceled. Let us 
now explain how to exploit this result even further. We start 
by defining the D x dD matrix L a ^ s p) = [(A s )U 1/2 ] aP . 
Clearly, the null space of this matrix is D(d — 1) dimen- 
sional. Let the Dd x D(d — 1) matrix V L with entries 
(as), 7 be a matrix of orthonormal basis vectors for this 
null space, which can be obtained from, e.g. the singular 
value decomposition of L, and thus satisfies LVl = and 
V[Vl = 1. We also introduce the notation V£ for the D x 
D(d — 1) matrix with components [Vj?] a , 7 = [Vj,]( a g), 7 . 
If we now group the (d — \)D 2 independent components 
of B in a D(d — 1) x D matrix x, we can use a parame- 
terisation B{x) given by B s (x) = l~ 1/2 V[xr~ 1 / 2 . One 
can check that this parameterisation satisfies the left gauge 
fixing constraint (l\E^} x ^ = since Vl contains only 

null vectors of L, and that B (x)GijB 1 (y) = \Z\ti[x^y], 
since the vectors in Vl are orthonormal. Up to the over- 
all diverging factor \Z\ that cancels in the LHS and RHS 
of Eq. (2), we have found a linear parameterisation B(x) 
for which the effective Gram matrix is the unit matrix. 
This same parameterisation cancels the last two terms in 
{dfip\H\^). The third term is still non-local, and requires 
the inversion of 1 — E. However, this is a pseudo-inverse 
as E has a single eigenvalue 1 and 1 — E is thus singu- 
lar. Let (K\ = (l\H£$(l - E)-\ We can safely replace 
(l\Hii by {l\Hii - h(l\, where h = (l\H^\r), since 
(Z|(l - E)- 1 = 0. Then, by replacing 1 - E with the 
non-singular matrix 1 — E + |r) (I \, we iteratively solve for 
the D x D matrix K from 

d 

K - ^(A s yKA s + tx[Kr]l = [{l\H%%\ - hi 

8=1 



with [{l\Hi£\ = Z stuv (st\h\uv)(A°Ayi(A"Av). 
Tracing this equation shows that trfi^r] = (K\r) = 
as required. Finally, we define the D(d — 1) x D tensor F 

d 

F=J2 (V / :) t / 1/2 C s V(A t ) t r- 1 / 2 

s,t=l 

d / d \ 

+ E(^) trl/2 J2^ n ° ts + KAS rl/2 ' 

s=l \t=l / 

where C st = ^2 uv (st\h\uv) A u A v . This defini- 
ton allows to write \\B t (x)\d^} — H\if;}\\ 2 = 
|Z|tr [x^x — x^F — F^x + constant] . This expression is 

minimized by choosing x = x* = F and thus A 1 = 
—\B{x*). Note that, thanks to the iterative solver, all steps 
can be performed in 0(D 3 ) computation time. 

Having now an explicit construction of A 1 , the simula- 
tion of time evolution with the TDVP now boils down to 
integrating a set of non-linear coupled differential equa- 
tions. The simplest numerical integrator is built on the Eu- 
ler method and proceeds as follows. 

1. Construct x* = F from the previous paragraph. 

2. Set A(t + dt) = A(t) - idtB(x*). 

3. Fix the gauge and norm of A by rescaling A. 

4. Calculate the energy and evaluate the step, change 
the time step dt if necessary. 

Step 3 is required since the gauge-fixing condition only 
fixes the norm and left eigenvector up to first order and 
higher order corrections are generally present. This sim- 
ple implementation is already useful for finding ground 
states through imaginary time evolution (dt — > —\dr). The 
TDVP then produces the best approximation to a gradient 
descent in the full Hilbert space which should be contrasted 
to a pure gradient-descent in parameter space such as in 
[11]. For real-time evolution, a simple first-order Euler in- 
tegrator does not inherit the symplectic properties of the 
differential equations and a more advanced integrator (see 
[8]) should be used. 

We now illustrate the power of our approach. Using 
imaginary time evolution with the simple Euler implemen- 
tation of the TDVP we've obtained a uMPS approxima- 
tion for the ground state of the S = 1 Heisenberg an- 
tiferromagnet. The TDVP stops when (d T ip\H\ip) = 0, 
which indeed signals a minimum in the energy expecta- 
tion value. Since the gradient has zero length at the min- 
imum, it automatically decreases in size as we approach 
it, and there is typically no need to reduce the size of the 
time step. This should be compared with the (i)TEBD case, 
where reduction of the time step, and thus automatic slow- 
ing down, is necessary to overcome the Trotter error. An 
ordinary laptop or pc allows one to find the ground state 
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TABLE I. First 24 Schmidt values of the D = 128 uMPS approx- 
imation for the ground state of the S — 1 Heisenberg antiferro- 
magnet. The degeneracy in the Schmidt spectrum as a result of 
SU(2) symmetry manifests itself, not by exploiting the symmetry, 
but rather by converging up to 'state tolerance' r\ = 10~ 10 . 



up to D = 1024 in less than one hour (without exploit- 
ing symmetries), resulting in a ground state energy den- 
sity e = -1.4014840389712(2) obtained with step size 
dt = 0.1. Since we can easily calculate the norm of the 
gradient as rj = \\x*\\, we can continue the evolution un- 
til ij has converged below a specified tolerance level. The 
convergence of the energy can be shown to be 0{rj 2 ) and 
can already be far beyond machine precision. This allows 
a much more accurate localization of the energy minimum 
than with the ordinary variational principle based on con- 
vergence of the energy, and is useful to e.g. obtain a very 
accurate convergence in the entanglement spectrum. The 
entanglement spectrum can offer valuable information but 
is not converged very accurately by other approaches (see 
[12] for an example). Table I shows how the first Schmidt 
values of the uMPS ground state for the Heisenberg chain 
at D = 128, which was converged up to rj = 10~ 10 , accu- 
rately reproduce the degeneracy according to half -integral 
spin representations. Note that we can also asses the er- 
ror of being confined to the manifold at any point in the 
evolution and derive from this a construction to optimally 
increase the bond dimension. Rather than starting from a 
random state at D = 1024, we can progressively build 
better approximations at larger D. Details are given in [8]. 

Using the time-reversal invariant numerical integrator 
discussed in [8], we can simulate a real-time evolution us- 
ing the TDVP equations. We start with the D = 128 uMPS 
ground state approximation of the XX-mode\ with mag- 
netic field \x = 1/2 along the z-axis, which is a critical 
model with non-zero magnetization {S z ) ^ 0, whereas 
(S x ) = (S v ) = due to the U(l) symmetry. We evolve 
this state according to the critical 5=1/2 Heisenberg 
antiferromagnet, so the expectation values (S x ' v ' z ) should 
be conserved due to the SU(2) symmetry. Comparative 
results for the TDVP implementation and a second order, 
translation-invariant TEBD implementation based on [13] 
are shown in Fig. 2 and illustrate that TDVP is much more 
capable of describing the evolution of conserved quantities. 

In this Letter we have introduced a new algorithm for 



simulating real and imaginary time evolution with (uni- 
form) matrix product states. The algorithm is shown to 
be globally optimal within the variational manifold, while 
conserving all symmetries in the system. 
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I This supplementary material discusses the following additional topics: 

• A review of the time-dependent variational principle (TDVP) and its properties by deriving it from an action 
principle. 

J>~j • Additional details about the implementation, including a description of convergence and error measures and a 

complete description of the time-reversal invariant numerical integrator. 

• A construction to dynamically increase the bond dimension of the uniform MPS. 

^vq • A generalization of the proposed strategy to the case of non-uniform MPS on finite lattices. 

^ REVIEW OF THE TIME-DEPENDENT VARIATIONAL PRINCIPLE 

In this section we review the general framework of the time-dependent variational principle, as can be found 
in [P. Kramer and M. Saraceno, Geometry of the Time- Dependent Variational Principle in Quantum Mechanics 
(Springer- Verlag, Berlin) (1981)]. Recall that the time-dependent Schrodinger equation (TDSE) can be derived from 
extremizing an action functional S{ip(t),il)(t)} — J* 2 dt L(vb(t), yb(t),t) with Lagrangian 

T3 i ! i . 

O 2 2 

O 



For the sake of brevity, we henceforth omit the time dependence of H and thus of L. Stationarity of the action under 
independent variations of \ip) and (ip\ in the full Hilbert space TL yields the TDSE and its complex conjugate. But 
when we only have access to a subspace or manifold M. C TL, we can still use the calculus of variations with respect 
to this action to define a time evolution of states \tp) e M; this is the essence of the TDVP. Hereto, we restrict to 
variations in the tangent plane of M at the point \vb). Assume that the manifold M can be parametrized as 

M = {|^(z)),zeC"}, (2) 

where we assume the dependence on the n complex parameters z % to be analytic and explicitly denote the anti-analytic 
dependence (^(z) | . Furthermore, we introduce the notation di for d/dz" 1 and c\ = djdz? , where we always use barred 
indices for the complex conjugate variables z. Requiring stationarity of S{z(t),z(t)} with respect to a variation 
z(i) — > z(t) + 6z(t) results in the following Eulcr-Lagrange equations 

iG ij (z(t),z(t))^(t) = <^(z(t))|ff|^(z(t))), (3) 

^ where G^j is the Gram matrix or overlap matrix of the tangent vectors of A4: 

G- lJ (z,z) = (ck^(z)\d^(z)). (4) 

We can also interpret this as the metric and — assuming linear independence of the tangent vectors — define the 
inverse metric as G^(z, z), such that G^(z, z)Gjfc(z, z) = 8 % k . The TDVP thus results in 

ii*(t) = G«(z(t),z(t)) (^{z{t))m{z(t))) (5) 

and its complex conjugate. In the main text, the Euler-Lagrange equations [Eq. (3)] were obtained geometrically, by 
looking for the coefficients z l (t) which minimize 

z\t)\d^{t)))-H\vb(z{t))) . 



2 



This minimization is obtained by the orthogonal projection of H \ip(z(t))) onto the tangential plane T^M, defined 

by 

T Z M = span{|<9iV(z)) ,i = l,...,n} . (6) 
The orthogonal projector Pt^m is indeed given by 

Pt z m(z,z) = \dnl>{z)) G«(z,z) (c^(z)| , (7) 

where the inverse of the Gram matrix appears in order to obtain P^m = Pt.m ■ 

Whereas Hamiltonian evolution in T~L is unitary and thus norm-preserving, this is no longer guaranteed for the 
projected evolution. In order to ensure norm preservation, we can define a modified Lagrangian L(ip(t),ip(t)) = 
L(i>(t),ip(t))/ (ip(t)\ilj(t)), which results in 

Z(z(t),z(i)) = i (fWdj-iFWdj) lniV(z(i),z(i)) H(z(t),z(t)) (8) 

where 

*<*.>-<*»*<.». ^-^^Mm 1 - 

Stationarity under variations z(t) + Sz(t) now results in the Euler-Lagrange equations 

iG,,(z(t),z(t))^(t) = ff,(z(i),z(i)), (9) 
and complex conjugates, where we have introduced the modified Gram matrix 

Gy(z,z) (^V(z)|^(z)>(^(z)|9^(z)) 



(z, z) = djdj In iV(z, z) 



iV(z,z) iV(z,z) 2 
(^(z)|a^(z)> (^(z)|V(z)) (V(z)|9,V(z)) 



(10) 



(V(z)|V(z)) (^(z)|^(z)) 2 
and the gradient of the normalized expectation value 

, _ (^(z)|g|^(z)) (^(z)|^(z)) _ (g#(z)|ff|y,(z)) (wwm*)) M*)\Hm*)) 

^(z,z)-d^(z,z)- ^.^ ff(z,z)_ W - Mz)) - (^(z)^)) 2 

These expressions can be easily interpreted. Under an infinitesimal variation, the norm or phase of a state changes 
if we move in the direction of \tp). Norm conservation is thus obtained when we subtract from every tangent vector 
\diip(z)} its component along \tp(z)) by replacing it with P (z,z) \diip{z)), where the projector P is given by 

. , mz)} (^(z)i 

P (z,z) = 1 - , (11) 

(^(z)|V(z)) 

and thus P (z,z) \diip(z)) = \ditp(z)) - \ip{z)) N&z)- 1 {^{z)\dii){z)) . We indeed find 

G y (z,z) = iV^z)- 1 (^(z)\P (z,z)\d^(z)) 

and 

H T (z,z) = Nfrz)- 1 (S^(z)|P (z,z)ff|V(z)) • 

If the manifold VW allows for norm and phase variations of states, i.e. if | V) € T 1 -^, then we can define the 
contravariant vector ip l such that ip l \ditp) — \ip). In this paragraph, we henceforth omit all arguments z and z for 
the sake of brevity. By definition we have that Pq \diip) ip l = and we can conclude that Gy has an eigenvalue zero, 
since Gyipj = = ip l Gjj, from which we immediately obtain the corresponding eigenvector. (Note that Gjj is a 
Hermitian matrix.) We now also define the covariant vector ifc = Gyip^ — (c\ip\ip) so that V'iV' 1 = (ip\ip) — N. With 
these definitions, we can write Gy = N~ 1 Gij — N^ 2 ^)^-. Even though Gy is not invertible, we can still define a 
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pseudo-inverse as G l1 = NG^ - , so that G^G- jk = 5\ - N~ VV>fc and GyG^ = S T k - N- 1 ^ . Since we can 

rewrite Hj as N~ 1 (S T k — N -1 ^-^ ) (dgip\H\ip), we are allowed to apply this pseudo-inverse to the Euler-Lagrange 
equations in order to obtain 

iz i (t) = G^(z(t),z(t))H J (z(t),z(t)). (12) 

In principle, the component of z(t) along the zero eigenspace of G can be chosen freely but, with the particular 
solution in the equation above, we satisfy ^^(t) = (ip\diip) z l (t) = 0, which is the required condition for norm (and 
phase) conservation. 

If the manifold Ai does not contain the freedom to change the norm and phase of a state, the modified metric G would 
have the same rank as the original metric G, which we assumed to be invertible. In particular, if \ip(z)) _L T Z A4, then 
Gij(z, z) = N(z, z)~ 1 Gij(z, z) and H T (z,z) = N(z, z) _1 (c\ipiz)\H\^p(z)) . The Euler-Lagrange equations following 
from S{z(t),z(t)} or from S{z(t), z(i)} are then identical. 

Finally, by defining for every pair of functions /(z, z) and g(z, z) a Poisson bracket 

{/• 9} = dif&djg - ru> J (13) 

we can write down the Euler-Lagrange equations as 

z l = i{H,z l } 7 f = \{H,T}. (14) 

For every operator O acting on W, we can define the expectation value 0(z, z) = (ip(z)\0\ip(z)) / (ip(z)\ip(z)) so that its 
time evolution is governed by O = i{H, O}. The manifold M. is thus a symplectic manifold. From the antisymmetry 
of the Poisson bracket we find {H, H} — 0, which implies that the energy of the state \ip) G M. is conserved under 
exact integration of the TDVP equations. 

The symplectic properties of the time-dependent variational principle also conserve other symmetries. Assume that 
the Hamiltonian is invariant under the action of a symmetry operator U (which should be a unitary operator), such 
that [H, U] = 0. In order to be able to transfer this symmetry to the manifold M. 1 we need to assume that for any 
state |V>( Z )) € M, the action of U is mapped to a new state |^(u(z))) = U |?A( Z )) <= -M.. Because of the unitarity of 
U, we have N(u(z), u(z)) = N(z, z), from which we obtain 

^#(z)G^ fc (u(z),u(z))9 ;U fe (z) - G- hl (z,z), (15) 

The condition [H, U] = also allows to conclude H (u(z), u(z)) = H(z, z) and thus 

^(z)ffj(u(z),u(z)) = H T (z,z), 
J ff J -(u(z),u(z))a i M J (z) = Hi(z,z). 

The metric and the gradient thus transform covariantly under the symmetry transformation and can be used to 
transform Eq. (9) into 

+i^#(z(t))G^ fe (u(z(t)),u(z(t)))^ U fe (z(t)) - as#(z(t))fl-j(u(S(t)),u(z(t))), 

and its complex conjugate. By using the injectivity of the map u(z), we can eliminate the Jacobians d^vP and d^v? in 
order to obtain the correct flow equations in terms of the new coordinates (u(t), u(t)). One case that is not covered 
by this general derivation is when U is an anti-linear operator, since u will then depend on z anti-holomorphically. 
Anti-linear transformations appear in quantum mechanics exclusively for time-reversal transformations. Let us denote 
R \ip(z)) — \ip(r(z))) with R the operator of an elementary time-reversal transformation. Because of the anti-unitarity 
of R and its commutation relation with the Hamiltonian (i.e. [H, R] = since we assume H to be time-reversal 
invariant), we obtain 7V(r(z)), r(z)) = N(z, z) = N(z,z) and H(f(z), r(z)) = H(z,z) = H(z,z), which yields 



and 



d l fHz)G %k (n^,r(z))d I r k (z) = %(z,z), (17) 

a i r?(z)H J (r(z),r(z))= J ff i (z,z), 
fl J -(r(z),r(z))a i r»'(z)= J ff T (z,z). 
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These relations convert Eq. (9) into 

+i^(z(i))G 5>fe (r(z(i)), r(z(t))Rr fe (z(t)) = H k (f(z(t)), r(z(t))Rr fe (z(i)), 
-ia 4 r fe (z(i))G Sj (r(z(t)),r(z(t)))^rnz(t)) = d# fe (zW)ff ¥ (r(z(i)), r(z(t))), 

or, by eliminating the Jacobian of the transformation, 

-iG w (r(z(t)2,r(S(t)))^H(z(t)) = ff,(r(z(i)), r(z(i))), 
+i^r?(z(t))C? 5ii (r(z(t)),r(z(t))) - ^(r(z(t)), r(z(t))) 

Note that the signs of the two equations have been switched, which is necessary to revert the time evolution of the 
new coordinates (r(i), r(t)). For a time-reversal invariant Hamiltonian H and a variational manifold M. that contains 
the time-reversed state R\ip(z)) G for each of its elements \ip(z)) <G A4, the flow equations of the time-dependent 
variational principle are also time-reversal invariant. This can be exploited in numerical integration schemes in order 
to construct symmetric schemes with improved stability (see next section). 

Returning to the case of linear symmetry transformations U, the expectation value U (z, z) is a constant of motion 
of the exact evolution according to the TDSE. For a general (discrete or continuous) symmetry, even when U \ip(z)) = 
|i/>(u(z))) e M, wc do not have that U\ip(z)) e T Z M. Consequently, U(z(t),z(t)) is not automatically a constant 
of motion of the TDVP evolution. But often U (z, z) does not represent an interesting quantity. However, when 
the symmetry operator U corresponds to a continuous symmetry generated by the Hcrmitian generator K, with 
[K, H] = 0, then the expectation value of K is an interesting constant of motion for the exact evolution according to 
the TDSE. We define a one-parameter family of transformations U(e) = exp(ie.K'). Since we require that for every 
state \ip(z)} in the manifold A4, U(e) \ip(z)) = \ip(u(z, e))) e A4, we can differentiate this defining relation with respect 
to e and set e = in order to learn 

\K |V(z)) - -^(z, 0) |9^(z)> = fc'(z) ■ (20) 

The action of K on a state |V>( Z )) thus has to be exactly captured in T Z M. We obtain for the part orthogonal to 
|V(z)) 

(dpl>{ti)\P (z,z)\d i il>(z))k\z) = i{d?l>(ti)\Po^^ 

iV(z, z)G^(z, z)F(z) = iiV(z, z)djK(z, z) 

and for the part parallel to \ip(z)) 

(^(z)|^(z))F(z) = i(^(z)|iv|^(z)) 
i>i(z, z)k i {z) = iN(z, z)K{z, z). 

The combination of these two equations results in: 

k\z) = iG 4j (z, z)djK(z, z) + iK(z, z)V J (z) 

If we now express H (z + ek(z), z + ek(z)) = (z, z) to first order in e, we obtain from the first-order expansion 

ie [illl(r 'i),l\ - diKG^djH + diHip'K - K^djH^j = 

(where we have omitted the arguments for the sake of simplicity) and thus, using \j} l diH = -ip J djH = H, we find 

{K 7 H} = 0, (22) 

so that any generator of a symmetry transformation of the Hamiltonian which can be exactly captured within the 
variational manifold produces a constant of motion K(z(t),z(t)) = K(z(0), z(0)). 

Finally, we can use the TDVP to simulate imaginary time evolution by setting t = —it. Imaginary time evolution 
according to the (imaginary) TDSE equation in the full Hilbert space % produces a gradient descent which converges 
to the true ground state |^o) for any initial state \ip) such that (V'olV') 7^ 0- There are no local minima in which the 
gradient descent can be trapped. Excited states correspond to saddle points. The imaginary TDVP evolution is the 



(19) 



(21) 
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best approximation to this gradient descent in "H; it is set apart from a gradient descent in parameter space (with 
Hi = di the gradient) by the appearance of the Gram matrix. Only when the Gram matrix G is the unit matrix, is the 
imaginary TDVP evolution equal to a gradient descent in parameter space. For the TDVP flow, the time derivative 
of the energy expectation value is given by 

^-H(z(t),z(t)) = -2d t HG^djH < 0. (23) 

The energy expectation value thus monotonically decreases until we reach the minimum, which is characterized by 
diH (z, z) = djH(z,z). Note that if we define the norm of the gradient as rj = [diHG 1 ^ d^H] 1 / 2 , then the rate of 
change of the energy expectation value is 0{rf). The energy thus converges quadratically faster than the state itself, 
which is a well-known result. 

Clearly, this whole section is applicable to the case of MPS and produces the results of the main text. In particular, 
the action of all symmetries for which the generator can be written as a sum of one-site terms can be exactly 
captured by the variational manifold -M u mps and thus produces constants of motion. The only difficulty in the MPS 
representation is its gauge invariance, which can be interpreted as an overparameterization and thus results in a set 
of tangent vectors which are not all linearly independent. But we have discussed in the main text how to overcome 
this difficulty. Since \diip) B l x is identically zero for all tangent vectors which consist of pure gauge transforms, any 
overlap or expectation value is also zero in the corresponding block. By using the linear representation B l (x), we are 
effective projecting all relevant quantities into their non-zero subspace, which is (d — 1)D 2 dimensional. 

DETAILS OF THE IMPLEMENTATION OF THE TDVP FOR UMPS 

In this section we discuss additional details of the implementation of the TDVP for uMPS. 

Choice of gauge 

In the main text we have discussed how to fix the gauge of the variations B l describing states in the tangent plane 
T^TVfuMPS- We call this gauge fixing constraint the left gauge-fixing condition, since it ensures that the left eigenvector 
and the largest eigenvalue (thus norm-preserving) of the transfer operator do not change to first order. Similar 
results could have been obtained with the following right gauge-fixing condition. We define the dD x D matrix 

iW = [r 1/2 {A s )^U (24) 

and define the (d — 1)D x dD matrix Vr that satisfies VrR — and VrV^ = 1. Thus, contains an orthonormal 
basis for the null space of Ftf. If we now also introduce [V^] 7iJ g = [Vr] 7i ( s /3) (with 7 = 1, . . . , (d — 1)1?, s = 1, . . . , d, 
and (3 — 1,...,D) then we find a different parameterization of the tangent plane as 

B s (x) = rWxVfr- 1 ' 2 , (25) 

where the independent parameters x have now been grouped inaDx (d— 1)D matrix. This parameterization satisfies 

the right gauge-fixing constraint E^ x ^\r) = and also produces B {x)GijB^ (y) = |Z|tr[a^y]. 

Up to this point, we have not elaborated on convenient gauges for A itself. A left orthonormalization gauge — 
where ^2 s =i(A s )^ A s = 1 anc ^ thus 1 = 1 — combines very well with the first-order conservation of the left eigenvector 
by variations B. Similarly, a right orthonormalization gauge — where r = 1 — combines very well with the right 
gauge fixing condition on B. However, both choices have also a numerical disadvantage. The representations B(x) and 
B(x) rely heavily on calculating l~ 1 / 2 and r -1 / 2 , which both appear in the different terms in the relevant expressions 
but never simultaneously (if worked out correctly). Hence, we would like to condition both matrices equally well at 
the same time. But if, e.g., I = 1, then r contains the eigenvalues of the density matrix of half of the chain, which 
are the square of Schmidt coefficients. If we try to obtain an accurate approximation where the smallest eigenvalues 
are of the order of machine precision, r is not determined accurately, it is ill-conditioned and the many operations 
with r -1 ' 2 can produce large numerical errors. Consequently, we have chosen to fix the gauge of A by requiring that 
I = r = A where A is the diagonal matrix containing the Schmidt coefficients. Starting from a general A, we first find 
the canonical form (A, T) [see G. Vidal, Phys. Rev. Lett. 98, 070201 (2007)] and then set A s = A 1 / 2 r s A 1 / 2 . This 
choice of gauge, which we call the symmetric gauge, evenly distributes the small eigenvalues of the density matrix over 
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the left and the right eigenvector of the transfer operator, resulting in a better conditioned algorithm. Nevertheless, 
these stability considerations are still the main reason that we cannot converge states up to machine precision (i.e. 
7] rts 1CT 15 ) at larger values of D where the spectrum of Schmidt coefficients contains many small values. 



Measures of convergence and error for uMPS 

Let's now reinterpret the TDVP evolution in order to assess both the convergence and the error of the TDVP 
equations. A measure of convergence is of course only useful for imaginary time evolution where we expect to end 
up in a steady state which is presumably the global minimum of the energy function H(A, A) and thus the best 
possible approximation of the ground state of H in the variational manifold .Mumps- Norm conservation is obtained 
by replacing the evolution vector H \if>) of the Schrodinger equation by 

P Q (A,A)H \i,(A)) = (i - \i,{A)) (1>(A)\)H \i>(A)) = (H - (ip(A)\H\ip(A))) |^) = (H - H(A, A)) , (26) 

where we have assumed that |^(j4)) is normalized to one. For a translation-invariant Hamiltonian H = X)«ez T n hT~ n , 
we introduce the notation (ij;(A)\H\ip(A)) = H(A, A) = \Z\h(A, A), where h(A, A) = (i/}(A)\h\ip(A)). Once again, we 
henceforth omit all arguments. If H is a nearest-neighbor Hamiltonian, so that h only acts nontrivially on sites zero 
and one, we have h = (l\H^\r) and we can write the exact evolution vector as 

(H - H) = T n v{ (• • • A s - 2 A s ~ 1 C s ° Sl A S2 ■ ■ ■ ) v R \. . . S ^ lSo s lS2 . . .) (27) 

with C st = J2t v =i ( s ' A^ 1 ~~ h\u, v) A U A V . The TDVP equation boils down to projecting this vector into the tangent 
plane TaM u mps, an d thus finding the set of coefficients x* such that 

B\x*)\drt) =J2 fn f E vl(---A s -iB s °(x*)A s +i---)v R \...S- lSo s +1 ...)\ (28) 

minimizes ||B l (a;) \diip) — (H — H) An algorithm for efficiently determining the set of the coefficients x* was 

explained in the main text. The additional term H \ip), not present in the main text, is non-essential, as it produces 
norm preservation which is automatically satisfied since (ip\diip) B t (x) = by construction. We now formally set 

B i {x*)\d^)=P TM {H-H)\^) 

as the vector that describes the TDVP flow. 

Let us now further elaborate on the properties of the TDVP flow for imaginary time evolution with uMPS. For 
imaginary time evolution, we can use the general result from the previous section to write 



±h(A(r),A(r)) = ±- *H(A(t),A(t)) = -±-\\Pr M (H H) |V>>|! 5 



^\\B^x*)\d^}\\ 2 = -tr[(xn\x*)]. (29) 



We obtain here some peculiar properties of quantum states of systems with infinite size (i.e. in the thermodynamic 
limit). Since the rate of change of the state (i.e. the norm of the gradient \\B l (x*) is proportional to IZI 1 / 2 , 

an infinite size uMPS always has zero overlap with the state it is converging to and with every other non-equivalent 
uMPS. This phenomenon is true in general for translation invariant states of a system of infinite size and was dubbed 
the (infrared) orthogonality catastrophe by Anderson [P. W. Anderson, Phys. Rev. Lett. 18, 1049 (1967)]. Indeed, 
let \ip(A)) and \ip(A)) be two normalized uMPSs (which requires that the largest eigenvalue of E\ and of E"~ is 1), 
where we also assume that both uMPS have a transfer matrix with a unique eigenvalue with magnitude one. We can 
then then compute the spectral radius p(E^) of E\, which is also know as the ground state fidelity per site d(A,A), 
since the total fidelity between the two quantum states scales as | (tp(A)\ip(A)) \ = F(A,A) = d(A, A)' Z L Then either 
d(A,A) < 1, such that (ip(A)\ip(A)) = 0, or d(A, A) = 1, in which case we can prove the existence of 4> € [0, 27r) 
and G e GL(D) such that A s — e 1<t 'GA s G^ 1 and \^(A)) ~ \%I){A)), provided that the boundary vectors «l and 
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are chosen identically. Two non- equivalent infinite size states (d(A, A) < 1) can however be locally similar and it is 
useful to define the local rate of change on the state as r] — ^/tr[(x*)tx*]. We then obtain for the rate of change 
of the energy density dh/dt — —i] 2 , which reproduces the quadratically faster convergence of the energy expectation 
value. Since convergence of the energy expectation value is the basis of all algorithms based on the (time-independent) 
variational principle, the corresponding state can still be far from the theoretically optimal state in the variational 
manifold. The same error propagates to all other expectation values of local operators. But in this imaginary time 
algorithm we can use the local change rj of the state as a convergence measure. At the point where the expectation 
value of the energy density is already at machine precision, we can further evolve the state in order to decrease rj to 
near machine precision and thus to very closely approach the theoretically optimal uMPS. Of course, the maximal 
accuracy that can be obtained depends on the conditioning of I and r, which are used in the construction of x* and 
B(x*), as discussed in the previous subsection. 

On the other hand, we can also assess the error between the theoretically optimal uMPS and the exact ground 
state. In fact, we can measure the error we make at any point in the TDVP evolution, both for real and imaginary 
time evolution. Note that the exact evolution vector is written in Eq. (27) in a form that is very similar to the tangent 
vectors, except for the fact that it acts non-trivially on two sites at the same time. By projecting onto the tangent 
plane, we are missing a part (1 — Ptm uM ps)(H ~ H) W) °f the exact evolution. We can thus assess the tendency of 
the exact evolution (H — H) \i[>) to take out of the manifold 7W u mps by calculating 



1 - Ptm uMFS )(H - H) |^)|| = ^AH 2 - B (x*)c\H - d i HB i {x*) + B {x*)G- l3 Bi {x*) 



AH 2 - B\x*) (drf)\drf) Bi{x*) ( 30 ) 
sJaH 2 - diHGvdjH. 



If this quantity is zero, then the evolution is exactly captured within our manifold, and the Hamiltonian can be 
thought of as effectively acting on a single site. In the final optimum where ||i? l (x*) \diip) |j = 0, the global first order 
error on the state is then given by the familiar expression \\(H — H) |i/>)|| = AH, which we now evaluate to find 



AH 2 = (t/j\(H - H) 2 \t/j) = {i>\T n {h-h)T m - n {h-h)T- m \i)) . 

n,m£Z 

With H being our nearest neighbor Hamiltonian and \ip) being our translationally invariant uMPS, we obtain 

AH 2 = |Z| ^ £ (M(h - h)f n (h - h)\1>) + 2(l\E% A (l E^E^r)^ . (31) 
The first term results in 

d 

(il>\(h-h)(h-h)\il>)= Yl (uMih-h^s^illA'A* ®A u A v \r), 

s.t,u.v—l 

d 

(i>\{h - h)f(h - h)\i/>) = v > w ^ h ~ h )f(h - tyf-^r, s, t) (l\A r A s A t ® A a A°A w \r), 

r,s,t,u,v,w—l 



for n — and n = 1 respectively, and in the complex conjugate of the last expression for n — — 1. Note that both 
quantities in the square root in Eq. 30 are proportional to |Z|. At any point in the evolution, we can thus define 
a local error on the state as e = |Z| _1 / 2 ||(1 — PtMumps)(H ~ ^0 IV')!!- I n imaginary time evolution, we converge 
towards the minimum where d l H = djH = 0. We thus find e = y/AH 2 /\Z\ in the optimal uMPS. The error m 
the expectation value of local operators is proportional to e, whereas the error in the expectation value of the energy 
density is proportional to e 2 . If at some point the error e exceeds a certain tolerance level we can expand the variational 
manifold by increasing the MPS bond dimension D. An optimal construction for this is described in the next section. 
This same construction also yields a better way to evaluate ||(1 — Ptm u mps)(H ~ H) IV')!!; which consists of a sum 
of quantities of comparable size which occasionally cancel each other. Large numerical errors can result from this 
cancellation. 
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Details of the numerical integration schemes 

When the original Hamiltonian is invariant under time reversal, i.e. RHRr 1 = H, then the TDVP equations 
have been proven to respect this time-reversal symmetry in the previous section, provided that R \ G M for every 
IV') € -A4. In quantum mechanics, the operator for time reversal is an antilinear operator. We now assume that the 
local basis is chosen so that R can simply be implemented as complex conjugation K, such that K\tp(A)) = \ij)(A)). 
This ensures that the time-reversed state of any uMPS is still a uMPS. Time-reversal invariance of the Hamiltonian 
then requires that all entries of h are real. Note that, in the standard choice of basis for a spin j system, R = e~ lvJy K, 
so that mere complex conjugation is in fact a combination of time reversal and rotation of the spin state over an angle 
7r around the y-axis. However, we ignore this subtlety and just refer to K as time reversal. 

In the case of imaginary time evolution, we can restrict to real representations for A, x* , and B(x*) whenever 
the Hamiltonian is time reversal invariant so that h has only real entries. A simple implementation based on the 
Euler method was sketched in the main text. Even though this is only a first-order method, and it introduces large 
(second-order) errors, this is not an issue for an imaginary time evolution due to the inherent stability of the approach. 
For imaginary time evolution with the TDVP we expect a monotonically decreasing energy expectation value and 
as long as the step size dr is small enough to reproduce this monotonic decrease, there is no need to change it. A 
decrease of time step should thus only be considered when higher-order effects result in an energy increase. The 
algorithm automatically slows down near the optimum since the norm Ic^V)!! approaches zero. The results 

for the S = 1 Hcisenberg antiferromagnet in the main text were obtained with a constant value dt = 0.1. This is in 
sharp contrast with TEBD implementations where there is a Trotter error associated with the time step dt, and the 
step size should be decreased as we approach the optimum. This decrease of step size produces an additional slowing 
down that has a significant effect on the total number of iterations required to converge the state. Finally, after every 
Trotter step, the resulting state is a matrix product state with increased bond dimension and we have to truncate it. 
In infinite-size systems, the only possibility is to use a local truncation based on the Schmidt coefficients, which is not 
globally optimal. All of these points have been illustrated in the main text. 

When using the TDVP to simulate real-time evolution, we are no longer able to restrict to a real representation. 
The quantities A, x* , and B(x*) become complex, even when the operator h has only real entries. Nevertheless, the 
algorithm sketched in the main text remains valid in principle. However, the large (second order) errors that are 
introduced by the Euler method can now accumulate in time. A prime indicator of this is a drifting expectation 
value (tp(t)\H\ip(t)) when the state \ip) is evolved in time according to a time-independent Hamiltonian H, whose 
expectation value should be conserved. The accumulation of systematic errors can be eliminated by implementing a 
numerical integrator for the Euler-Lagrange equations that respects the symplectic structure of the TDVP. However, 
the structure of the TDVP is much more complicated that the typical structure of classical dynamics with a separable 
Hamiltonian H(q,p) = T(p) + V(q). In particular, the relation H(A,A) is highly nonlinear and not separable. 
Consequently, none of the typical symplectic algorithms from classical dynamics can be applied to the TDVP. 

If a set of differential equations is invariant under time reversal, it is a good policy to devise a numerical integration 
scheme that respects this time reversal symmetry. Numerical integration schemes that respect time reversal symmetry 
are called symmetric and share many nice properties with symplectic integration schemes such as a stable long-time 
behavior, a linear growth of the global error and a near-preservation of first integrals [E. Hairer, C. Lubich and 
G. Wanner, "Geometric Numerical Integration: Structure preserving algorithms for ordinary differential equations", 
Springer Series in Computation Mathematics 31, Springer (2002).]. The following paragraphs describes the details 
of a second order numerical integration scheme that respects time-reversal symmetry, although it can of course also 
be applied to Hamiltonians which are not invariant under time reversal, in which case it is a simple second-order 
numerical integrator. At any point, the calculation of the projection of (H — H) \tp) into the tangent plane TaM in 
the point \i/)(A)) (i.e. the determination of x*) can be performed by the algorithm that was outlined in the main text. 

The main problem one encounters when trying to construct a symmetric integrator for differential equations on 
manifolds is that the tangent plane at \ip(A(t))) and \ip(A(t + dt))) are different. Most algorithms thus start with 
the determination of a midpoint \tl>(A)) such that the projection of Pt a m IVK^W)) + Pt a m l^i^t + dt))) = 0, with 
TjAi the tangent plane in the midpoint. This relation only specifies the location of the midpoint, but not the relation 
between A(t) and A(t + dt). In principle, the midpoint can be combined with any integration scheme, where we 
use the midpoint to calculate the variation B(x*). If we thus use the simple Euler step, we obtain the additional 

relation Pt a m \ip(A(t + dt))) - Pt x m \ip{A(t))) =idtP TAM (H - H) \ip(A)) where H = H(A, A) is the energy in the 
midpoint. This is_ sketched in Figure I. This specific combination of midpoint with Euler step immediately tells_us 
that A(t + dt) = A + idt/2B(x*), where x* is the set of parameters such that B' L (x*) \diip(A)) = P TaM (H-H) \ip(A)). 
Note that the representation of the variations B itself depends on the current value of A, a subtlety which we capture 
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FIG. 1. Sketch of the location of the midpoint to be used in a symmetric integration scheme on a manifold 



by simply introducing a suitable notation for the argument x. A similar reasoning leads to the conclusion that the 
midpoint is implicitly defined by A(t) = A — i dt/2 B(x*). However, we also have to take into account the normalization 
preservation and gauge fixing that is applied in every step. Instead of setting A(t + dt) = A + idt/2 B(x*), we find a 
scalar constant c and a D x D matrix G £ GL(D) such that A s (t + dt) = cG[A s + idt/2B s (x*)]G~ 1 , where c and G are 
chosen such that has 1 as largest eigenvalue, and satisfied a prescribed gauge fixing condition, such as the 

symmetric gauge defined in the previous subsection. Analogously, we also have to look for (different) c and G in the 
defining relation c~ 1 G~ 1 A s (t)G = A s — i dt/2 B s (x*), which are chosen such that c~ 1 G~ 1 A s (t)G — A s is compatible 
with the gauge fixing constraints that are built into the representation B(x). Since every B(x) satisfies (l\E^ x ^ = 0, 
we have to tune c and G such that 

(T\Ef G ~ lA{t)G - A = o => (Tg- 1 ^ = c(lG _1 |. 

Put differently, c is the largest eigenvalue of E^® and G is chosen such that (IG~ 1 \ is the corresponding left eigenvector, 

where (l\ is the eigenvector of E~ corresponding to eigenvalue 1. Since we cannot solve the resulting implicit relation 

c~ 1 G~ 1 A s (t)G = A s — idt/2 B s (x*) exactly, we have to devise a numerical scheme to determine A. The resulting 
algorithm doesn't satisfy time reversal symmetry exactly, but only up to the accuracy of the numerical determination 
of A, which can be near machine precision. 

We can try to solve the implicit relation for the midpoint by a simple error correct strategy. We choose as an initial 
guess A <~ A(t) + idt/2 B(x* (t)), where the 'similarity sign' is used to indicate that A has already been transformed 
in order to satisfy norm and gauge fixing constraints. Having a guess A ni we can iteratively try to improve it as 
follows. We calculate the difference dA n = c~ 1 G~ 1 A s (t)G n — A s n + \dt/2B(x* n ), where c„ and G n are chosen such 
that (InG^lE^^ — c n {l n G~ 1 \. We then set A n+1 ~ A n + dA n and repeat this process. At any point in the iteration, 

we can measure the size of the correction as \\dA n \diip(A n ))\\ = \ l E\ 1 / 2 (l n \E'^ n \r n ) 1 ^ 2 . As argued in the previous 

subsection, we can safely omit the overal \ r L\ 1 ^ 2 in order to obtain a local measure £ = {l n \E'^ n \r n ) 1 /' 2 . When £ dives 
below a tolerance level that can be chosen near-machine precision, we can stop the iteration. When the chosen time 
step is not too big — for example dt w 0.01 — this algorithm converges in a few (less then 20) iteration steps. Better 
strategies in terms of higher order iterative solvers for non-linear equations can be devised. 
The general outline of an algorithm for real time evolution is thus: 

1. iteratively determine the midpoint from A(t) ~ A — i dt/2 B(x*) 

2. set A(t + dt) - A + idt/2B(x*) 

Note that all operations can be implemented in 0(D 3 ) computation time. We of course use iterative eigensolvers to 
determine the eigenvalues and eigenvectors of i?-^, £/~ and E^^}. 

The midpoints \A) can in fact be interpreted as \A(t + dt/2)). Thus the algorithm produces twice the resolution as 
initially requested. However, it is only (approximately) time reversal invariant after an integral number of steps dt. 
Furthermore, since the (backwards) Eider method is used to step from A(t) to A(t + dt/2), the error in this step is 
expected to be 0(dt 2 /4). Similarly, the error in the step A(t + dt/2) to A(t + dt) is expected to be of the same order. 
Nevertheless, the resulting step from A(t) to A(t + dt) is correct up to second order, and the error is actually 0(dt 4 ) 
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because odd-powered effects are forbidden by the symmetry of the construction. Higher order errors are obtainable 
by combining the midpoint construction with more advanced Runge-Kutta schemes. 

We can once again compare this implementation of the TDVP with a real time iTEBD simulation. Not only 
does the TDVP-based algorithm for real time evolution have the same advantages as the imaginary time algorithm 
- conservation of translational invariance and internal symmetries — the (approximate) time-reversal invariance 
makes the algorithm extremely stable over longer simulation times. In principle, the first step in a real time iTEBD 
implementation is even better, because evolving over a small time using a Trotter decomposition is a symplcctic 
operation. However, since iTEBD takes the state outside the manifold of MPSs with fixed bond dimensions, this step 
is followed by a truncation that breaks both the symplectic symmetry and the time reversal symmetry. 

After having applied the real time evolution algorithm to evolve a state \ip(A)) at t — to \ip(A(t))) at time t for a 
total time tf, we can explicitly apply the time reversal operator K to the final state \ip(A(tf))} and use the resulting 
state \tp(A(tf))) as the initial value for a new time evolution over a total time tf, resulting in states \tp(A'(t)} at time 
t with thus A'(0) = Ajtfj. We should then compare the states K \ip(A'(t f - t))) = \ip(A{t))} [with A(t) = A'(t f -t)] 
to the states \i/j(A(t))). By using the results from the previous subsection, we know that the equality between two 
uMPS \ip(A)) and \ip(A)) can be measured by computing the ground state fidelity per site d(A, A) = p(E^), which 
is one for equivalent uMPS and will be smaller than one for states that differ. Fig. 2 compares the time-reversal 
symmetric integrator in this subsection with results a TEBD-based algorithm, for the matrices A(t) were taken from 
the results in Fig. 2 in the main text. 

Klxln -, 1 
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FIG. 2. Comparison of time reversal invariance in a TDVP (dashed lines) and TEBD (dotted lines) simulation. Illustrated is 
(one minus) the ground state fidelity d(A(t), A(t)), where A(t) = A'(t f - t) (t f = 10, dt = 5 x 1CT 3 , C = 10~ 10 (see text)) 



DYNAMICAL EXPANSION OF THE VARIATIONAL MANIFOLD 



Both for real and imaginary time evolution, we have introduced the local error measure e = 



-1/2 



11(1 - 



Ptm u mps)(H ~ H) WW that captures the tendency of the exact evolution to move away from the manifold .M u mps- 
If this quantity exceeds a given tolerance value, we might try to reduce the error by expanding the variational class. 
For uMPS with bond dimension D, we can expand M u mps(d) by increasing the bond dimension to some value D. 
If the state \ip) € M u mps(d) at some point is a uMPS with D x D matrices A s , we can try to better approximate 
the exact evolution by embedding at the next iteration this state into a larger manifold -M uMPS ^ by defining new 

D x D matrices A s with D > D and 



A' 




dA s 



where the variation dA s should be proportional to the chosen time step dt, and is given by 

dA s = 



dAf) dAf)! 



The first order effect of this evolution step is given by 

d 



|#o) = ^T"[ J2 vt(---A s -i(dA 00 y°A s +i---)v R \...s- lSo s +1 ...) 



(32) 



(33) 



(34) 



nGZ 



which is equivalent to the evolution in the original manifold. The newly introduced degrees of freedom do not seem 
to appear at first order, which is a consequence of the intrinsic nature of the TDVP not to leave the manifold in 
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which it is defined. However, by choosing dA m and dA w proportional to dt 1 ^ 2 , they do in fact generate a first order 
contribution 




(• • • A a -*A'- 1 (dA 01 ) a °(dA 10 )'+iA a +* ■ ■ ■ )v R |. . . s _ 2S _ lSo s+i S+2 • • •) 



(35) 



Clearly the increased bond dimension allows for an action on two neighboring sites, which increases the entanglement 
entropy from some initial value smaller than log(Z?) to some new value that is smaller than the new bound log(D). 
Note that dAf x can only appear in first order if it is of order O(dt ) = 0(1). At this point, dA^ can appear any 
number of times as a string in between dA^ and dA s w . Thus dA^ induces a non-trivial affect on any number of 
sites, with a minimum of three sites. It is thus not useful for nearest-neighbor Hamiltonians, and is henceforth set to 
zero: cL4 n = 0. We show below that, by using the correct gauge conditions, the overlap between the exact evolution 
vector for a nearest neighbor Hamiltonian and such a vector with action on more than two sites is exactly zero. The 
freedom in dA s lx can possibly be useful in Hamiltonians with long range interactions. 

If we try to generalize the geometric strategy of the TDVP, we should look for optimal matrices eL4g , dA^ and 
gL4| such that 



||di(tf-iO|V)-|#o>-MV>i 



(36) 



is minimized. This is a complicated expression that couples the three sets of parameters cL4g , dA^ and cL4f - Here 
too, we can apply infinitesimal gauge transformations 



G 



1 + eXoo e 1/2 Xoi 
e 1 ' 2 X w 1 + eXu 



from which we obtain up to 0(e) 



G 



~A S 


0" 


G- x = 


~A S 


0" 















+ dA x with dA x 



e[X 00 , A s ] + f {X 01 X W , A s } -e 1 / 2 A s X l 



01 



z l l 2 X w A s 



-eX w A s X { 



01 



Any dA is thus gauge equivalent to dA + dAx- Note the absence of Xn at order e. It is also important to note 
that exploiting gauge freedom to reduce the number of independent parameters mixes the off-diagonal blocks with 
the diagonal blocks. Consequently, we cannot generally treat the optimization with respect to dA 00 independently 
from the optimization with respect to cL4 i and dA w . For evaluating (dipo\dipo) and (dipo\(H — H)\ip), we can use 
the results from the main text. The corresponding relations involving \dtpi) are given by 



+ (l\E^(l E)-'E d //\r) + (l\E d / A Hl E)~'EX\r) + |Z|(i|^|r)(/|^|0 



(37) 



(dip \dipi) 



(l\E 



AdAoo 1 



\r) + mZlAV) - 2{l\E^ A Jr){l\E^\r) 



+ ('!<,(! - E)-*E%£\r) + (l\E d //(l - Ey'EXjr) + |Z|(/|£7^ 00 |r)(i|^|r) (38) 



and 

= |Z| 



{l\E c dAl \r) + {l\Ej% lA \r) + {l\E c J Al \r) 

+ (l\E^(l E)^E c AA \r) + mZAl E)-^EXV) 



(39) 



where dAf = dA^dA^. While finding the simultaneous optimum for dA o, dA i and dA w from these equations 
appears to be a very complicated task, we can exploit the gauge freedom to eliminate most terms. We can impose 
(l\E A Ao1 = and E A Aw \r) = by exploiting the D x (D - D) parameters in X 01 and the (D - D) x D parameters 
in Aio respectively. This automatically removes all non-local terms from (di/^i M^i) an d {dtpi\H — H\ip), i.e. 



(diP 1 \diP 1 } = \Z\(l\E d d il\r), 



{d*h\H-H\ 



E C dAl \r). 
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In addition, this choice of gauge also renders (ip\dipi) — and (dip \dipi) = for any dA no . This choice of gauge thus 
ensures that \di(>i) is orthogonal to both the state \ip) as well as to the complete tangent plane TaM u mvs(d) an d omv 
captures effects that act on two sites at once. These effects cannot be captured by \dipo) € TaA4 u mps(d)> where only 
a single site is modified. The optimization problem for cL4oo completely decouples from the optimization problem for 
dAoi and dA w , and can be solved as described in the main text. Note that we still have the D 2 — 1 parameters in X no 
together with global norm conservation (tp\dipo) + (V'MV'i) = (V'MV'o) = to fix the gauge of dA 00 as described in the 
main discussion. Furthermore, as stated above, with this choice of gauge for dA 01 and dA 10 , the tangent vectors with 
nonzero choice of dA\i are orthogonal to the evolution vector (H — H) for any nearest neighbour Hamiltonian. 

For the optimization problem in dAoi and cL4io, we can use a parameterization that is based on the definitions in 
the first subsection of this section. We define Bq 1 (x) = l^^V^x with x a (D — l)d x (D — D) matrix of independent 
components, and analogously Bf (y) — yV^r -1 ' 2 with y a (D — D) x (D — l)d matrix of independent components. 
We then define Bf t (x,y) — BQ 1 (x)B{ (y) in order to find 

m B B %^\r) = tr [(xy)(xy?] , (!|< (l , v) W = tr [G{xy?] , 

with G = Ef. t= i(^) t;1/2 C sirl/2 (^) t - Since we have to minimize 

tr [(xy)(xy)i] - tr [G(xy?] - tr [(xy)&] , 

we are looking for the optimal matrix xy of rank D — D such that \\xy — G\\ is minimized, with ||-|| the Hilbcrt- 
Schmidt norm. The best approximation (x,y)can be found by performing a singular value decomposition of the 
(d — 1)D x (d — 1)D matrix G and retaining the largest D — D singular values. If D = dD, we can always capture 
one step in the action of the Hamiltonian exactly, which is also the case with the TEBD algorithm. 

This construction is very useful in combination with imaginary time evolution in order to correctly converge a 
random initial state to the optimal state within the variational manifold. Rather then starting from a random uMPS, 
it is better to use a previous optimal uMPS A at a smaller bond dimension D and_combine it with the construction 
above to form an initial state A for the simulation at the new bond dimension D. Note that, if A represents the 
optimal solution at D, we can choose dAf )0 = 0, or thus \dipo) = 0. This approach avoids a particular problem of the 
imaginary time evolution with the TDVP: For large D, the current TDVP implementation is susceptible to converging 
some random initial states to a state where E\ has more then a single eigenvalue 1. At this point, the formulas for the 
gradient and the gram matrix in the main text are no longer valid and incorrect application of these formulas often 
results in the algorithm getting trapped. While we could of course implement more advanced formula's that take this 
scenario into account and with which the implementation would be able to continue the convergence process, using 
the aforementioned construction avoids the occurrence of this problem completely. In addition, very few iterations 
are needed at the large values of D, since the error e of the previous result at slightly smaller D is already very small. 
Clearly, this is the preferred approach. 

In principle, the same construction can be used in combination with real time evolution. After quantum quenches, 
the entanglement entropy typically increases and so does the error measure e. When e increases beyond a certain 
tolerance level, we can choose to use the above construction to dynamically increase D, which brings e back to an 
acceptable level. However, since this whole step is only well defined at first order, it cannot be included in a higher 
order numerical integrator. The dynamic expansion of the variational manifold should thus be used carefully in real 
time evolution, since it also breaks the time reversal symmetry. 

Finally, we can also use the results from this section for another purpose. Since \dip ) captures the projection of 
(H — H) \ip) in the tangent plane TA4, it is indeed equal to B l (x*) \diip) = Pt a m uM ps(o) {H ~ H) l^)- ^ we choose 
D = dD, then \dipi) can exactly capture the remnant (1 — Pt a m uM ps(o))(H — H) \tp) that is orthogonal to the tangent 
plane TaM u mps(d)- We can use this to solve the following problems. If the exact evolution is almost exactly captured 
in the tangent plane TaMuMPS(d), the error measure e has to be calculated as the square root of the difference of two 
'equally large' numbers AH 2 and IcW'^))!! 2 - It is then numerically much more convenient to calculate e as 

e = |Z| _1 / 2 || I^V^)!! = m^dAl l 7 ") 1 ^ 2 ' which can completely be constructed without having to subtract large numbers 
with small differences. In addition, the expression for AH 2 given in Eq. (31) also contains 4 contributions which 
might almost cancel each other if \ip) is very close to an cigenstate of H. This too can produce large numerical errors, 
and it is better to calculate AH as AH = \/ (V'olV'o) + (V'iIV'i) where the square root contains two positive numbers. 



13 



TIME-DEPENDENT VARIATIONAL PRINCIPLE FOR GENERIC MPS 



We conclude this supplementary material by investigating whether new efficient algorithms can be obtained from 
applying the TDVP to generic (non-uniform and finite-size) MPS. The TDVP equation for finite size MPS was already 
written down in a different format in [J.J. Dorando, J. Hachmann, G.K.L. Chan, J. Chem. Phys. 130, 184111 (2009)], 
but the resulting equation was not further investigated. We define a MPS with open boundary conditions as 

d 

\i>[A{n)])= J2 vlA s Hl)A^(2)---A^(L)v R \s lS2 ... SL ) (40) 

with site-dependent matrices A s (n) having site-dependent dimensions D n _\ x D n . It is always possible to absorb 
into A Sl (l) and set Dq = 1; similarly we can absorb into A SL (L) and set Dl — 1- This completely eliminates the 
need for the boundary vectors, and they are henceforth omitted. We can now define the series l(n) and r(n) (with 
n — 0, . . . , L) of site-dependent D n x D n density matrices for the auxiliary system through 1(0) = 1 = r(L) and 

d d 
l{n) = ^{A s {n))n{n-l)A s {n), r(n) = ^ A s (n + l)r(n + l)(A s (n + l))t. (41) 

s=l s = l 

The norm of the state is then given by r(0) = l(L) = (l(n — l)|£^"j|r(n)), which we require to be one. 

If we want to apply the TDVP to the class Mmps — {\4>[A(n)])} of MPS of length L with fixed bond dimensions 
{D ni n = 1, . . . , L}, we have to vary with respect of all entries in every tensor A(n). We now denote a general variation 
of A as B l with a collective index i — (a, s, 0, n), such that the tangent space T^^AImps is spanned by the states 

L d 

B i \d i ^) = J2 E A s >(l)---B s "(n)---A s -(L)\ Sl ... Sn ... SL ). (42) 

n=l {s n } = l 

We can easily calculate the metric Gy = (drip\djij>} , but it is a complicated expression coupling all variations, so 
it would seem like an impossible task to invert this LD 2 d x LD 2 d matrix. However, as for the case of uMPSs, we 
can exploit the gauge invariance to remove all non-local couplings, i.e. all contributions containing variations at 
different sites. Combining the particular choice of gauge fixing conditions which has this effect with a well chosen 
linear parametrization yields an effective metric that is the unit matrix. 

Firstly, we note that by choosing B s (n) — A s (n)/L 7 we have B % \diip(A)) = \^p{A)). We thus have to project onto 
norm-preserving variations. Secondly, we now have a site dependent gauge freedom A s (n) <— G(n — l)A s (n)G(n) _1 
with G(n) e GL(D n ). From the first order effect of an infinitesimal transformation G(n) = 1 + eX(n), which is 
A s (n) <- A s (n) + eB s x {n) with B s x (n) = X(n - l)A s {n) - A s {n)X(n), we conclude that B x \diip(A)) = as can 
easily be checked explicitly. Thus any variation B l is gauge equivalent to B l + B l x for all {X(n) e C D " xD "}. 

We can fix the gauge of the variations (and ensure norm conservation) by requiring that cither (Z(n — 1)\E^™} = 
or E^2)\ r ( n )) = f° r cver y n , which we call left and right gauge fixing conditions. Both choices reduce the metric to 



B'~ l G-^ = J2(l(n l)\E^ } \r(n)). (43) 



If J? is a nearest neighbor Hamiltonian H = X)n=i M n ) with h(n) acting non-trivially only on sites n and n+ 1, then 
the corresponding right hand side of the Euler-Lagrange equations reduces for the case of left gauge fixing to 



B T (tk\H-H\i,) = J2 



n=l 



6(n<L){l(n-l)\E C B ^ )A{n+l) \r(n + l)) 

+ 9{n > W(n 2)\E^:l] B{n) \r(n)) (44) 

n-2 / n-1 \ i 

+ 0{n > 2) ^ (Km - l)\E C A $ A{m+1) ft E$f> E A B ^\r(n)) 



m=l \k=m+2 



14 



where 9 is a discrete Heaviside function that yields one if its argument is true and zero otherwise, and where C st (n) = 
£*,„=! ( s > - &(»)!«>«> A"(n)A"(n + 1) and 

d 

h(n)= Yj (s,t\h(n)\u,v)ti[l(n- l)A u (n)A v (n + l)r(n + l)(^ s (n)A*(n + 1))*] . (45) 

In case of right gauge fixing, the last term in Eq. (44)would be replaced by a term containing all contributions of the 
Hamiltonian acting to the right of B(n). These terms are familiar from the variational sweeping algorithm for MPS, 
and it is well known how to construct them efficiently and iteratively. Clearly now, all variations B(n) decouple from 
variations B(n') at different sites n' ^ n. 

Analogously to the case of uMPS, it is easy to find a parameterization of B[x] depending on L matrices x(n) 
(n = 1, . . . , L) of size (cLD„_i — D n ) x D n , such that B[x](n) depends only on x(n) and such that the effective Gram 
matrix becomes the unit matrix. We first define the D n x <il?„_i matrix 

[i(n)]a,(.« - [A'(n)H(n - l) 1/2 ] a , fi (46) 

and then construct a dD n -\ x (dD n _\ — D n ) matrix Vl(ji) that contains an orthonormal basis for the null space of 
L(n), i.e. L(n)Vt(n) = and Vl(7i)^Vl(™) = 1, for all n = 1, . . . , L. We then use the representation 

B s [x](n) = l(n - l)- 1/2 V[x{n)r{n)- 1/2 (47) 

in order to obtain 

L 

B^x^B^y] = tr [x(n^y(n)] . (48) 

n=l 

Let us now outline all the necessary steps in a single iteration of a simple Eulcr implementation for imaginary time 
evolution: 

1. Define K(0) = (scalar) and K(l) = (D\ x D\ matrix) and compute the D n x D n matrix K(n) as 

d d 

K(n) = ^ A*(n)U s (n - l)H(n - 2)C st (n - 1) + ^ A s (n)^K(n - l)A s (n) 

s,t=l s=l 

for n = 2,...,L-l. 

2. Define the (dD„_i — D n ) x D n matrices F(n) as 

d 

F(n) =6{n < L) ^ V£(n)tj(n - lJ^C^nMn + l)^*(n + l)V(n)- 1 / 2 

s,t=l 
d 

+ 0(n > 1) X! F i _ l) _1/2 A*(n - l)t/( n - 2)C ts (n - l)r(n) 1 / 2 (49) 

s,t=l 
d 

+ 0(n > 2) ^ VZ(n)ti(n - \y x l 2 K (n - ljA'^rW 1 / 2 . 

3. Set z*(n) = F(n). 

4. Take a step A s (n,i + di) = ,4 s (n) - d<B s [x*](n). 

5. Set 1(0) = 1 and compute l(n) = ^2^ = i{A s (n))U(n — l)A s (n) for n = 1, . . . , L. Renormalize if necessary. Set 
r(L) = 1 and compute r(n) — £ s=1 A s (n + l)r(n+ l)(A s (n + 1))* for n = L — 1, . . . , 0. Also compute h(n) and 
change gauge fixing if necessary. 

6. Compute H — J2n=i M n ) an d evaluate the step. Change dt if necessary. 
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Every step can be computed with computation complexity 0(LD 3 ), just as in the standard sweeping algorithm. 
However, unlike in the sweeping algorithm, the determination of the change A s (n) <— A s (n) — dtB s [x*](n) — even 
though only valid for dt not too large — is globally optimal and thus includes all correlations from the changes at 
every other site in the lattice. In addition, the variation B s [x*](n) does not require to (iteratively) solve an eigenvalue 
equation but is determined directly from matrix multiplications of 0(D 3 ). Since the iterative solution of the the 
eigenvalue problem is the dominating contribution to the computation time for large D in the traditional sweeping 
algorithm, it is interesting to see whether this implementation can outperform the traditional sweeping algorithm in 
convergence speed. Clearly, all the other aspects of this supplementary material (a time-reversal invariant integrator 
for real time evolution, dynamical expansion of the variational manifold, error measures, . . . ) can also be applied 
to the non-uniform, finite-size lattice. Note that, for lattices with periodic boundary conditions, the gauge freedom 
cannot be exploited to fully decouple variations at different sites, resulting in a much more complicated algorithm, 
which is also the case with other DMRG-related approaches for lattices with periodic boundary conditions. 



